#include "mvn.h"
#include "MersenneTwister.h"
#include <fstream>
#include <sstream>
#include <iostream>

//Global variables
size_t size;		//Size of the matrix and vectors
size_t Nmax;		//Max times of iterations
Matrix<double> *pInput(NULL);				//Input convariance matrix
vector<double> *pVecA(NULL), *pVecB(NULL);	//Input vectors
double eps;
double alpha(4.5);

bool ReadInput(string);						//Read file for initialization
void Release();								//Release memories
 
int main(int argc, char *argv[])
{
	if(argc != 2)
	{
		cout<<"Usage:[filename] [InputFileName]."<<endl;
		return -1;
	}
	//Read datas from file for initialization
	if(!ReadInput(argv[1]))
	{
		cout<<"Input error!"<<endl;
		return -1;
	}

	Matrix<double> C = CholeskyFactor(*pInput);			//Calculate Cholesky Factor C of the input covariance matrix
	
	//Intsum: final result
	//Varsum: temp variables
	//N: number of iterations
	double Intsum = 0.0, Varsum = 0;
	size_t N = 0;

	//Three vectors: d,e and f, used in the algorithm
	vector<double> d(size),e(size),f(size);
	d[0] = Phy((*pVecA)[0]/C[0][0]);
	e[0] = Phy((*pVecB)[0]/C[0][0]);
	f[0] = e[0] - d[0];
	
	//Random number generator
	MTRand random;

	//Start the algorithm
	//The main loop
	do
	{
		//Generate random numbers w1,w2,..., size-1 in total
		vector<double> w(size-1);
		for(vector<double>::size_type i=0; i<w.size(); ++i)
		{
			w[i] = random.rand();
		}

		//for 2,...,size
		vector<double> y(size-1);
		vector<double> tmpSum(size);
		for(size_t i=1; i<size; ++i)
		{
			y[i-1] = InvPhy(d[i-1] + w[i-1]*(e[i-1]-d[i-1]));

			tmpSum[i] = 0.0;
			for(size_t j=0; j<i; ++j)
			{
				tmpSum[i] += C[i][j]*y[j];
			}
			//Update vectors d, e and f
			d[i] = Phy(((*pVecA)[i]-tmpSum[i])/C[i][i]);
			e[i] = Phy(((*pVecB)[i]-tmpSum[i])/C[i][i]);
			f[i] = (e[i]-d[i])*f[i-1];
		}

		++N;
		double sigma = (f[size-1]-Intsum)/N;
		Intsum += sigma;
		Varsum = (N-2)*Varsum/N+sigma*sigma;
	}
	while(N != Nmax);

	cout<<"Result:\t"<<Intsum<<endl
		<<"Iteration time:\t"<<N<<endl;
	
	//Release memories before exit
	Release();

	return 0;
}
	
bool ReadInput(string fileName)
{
	//Read the input file
	ifstream file(fileName.c_str());
	if(!file)
	{
		cout<<"File '"<<fileName<<"' not fould!"<<endl;
		return false;
	}

	string line;
	//Read lines one by one
	while(getline(file,line))
	{
		//'#N' means the size of the matrix and vectors
		if(line == "#N")
		{
			getline(file,line);
			istringstream iss(line);
			iss>>size;
			pInput = new Matrix<double>(size);
			pVecA = new vector<double>(size);
			pVecB = new vector<double>(size);
		}
		//'#MAXPTS' means the max number of iterations
		else if(line == "#MAXPTS")
		{
			getline(file,line);
			istringstream iss(line);
			iss>>Nmax;
		}
		//'#ABSEPS" for error check, not used
		else if(line == "#ABSEPS")
		{
			getline(file,line);
			istringstream iss(line);
			iss>>eps;
		}
		//'#A': Datas for input vector A
		else if(line == "#A")
		{
			getline(file,line);
			istringstream iss(line);
			for(size_t i=0; i<size; ++i)
			{
				iss>>(*pVecA)[i];
			}
		}
		//'#B': Datas for input vector B
		else if(line == "#B")
		{
			getline(file,line);
			istringstream iss(line);
			for(size_t i=0; i<size; ++i)
			{
				iss>>(*pVecB)[i];
			}
		}
		//'#Cov': Datas for input convariance matrix
		else if(line == "#Cov")
		{
			for(size_t i=0; i<size; ++i)
			{
				getline(file,line);
				istringstream iss(line);
				for(size_t j=0; j<size; ++j)
				{
					iss>>(*pInput)[i][j];
				}
			}
		}
		else
		{
			getline(file,line);
		}
	}

	return true;
}

void Release()
{
	delete pInput;
	delete pVecA;
	delete pVecB;
}